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1  Objectives 

We  seek  to  produce  fast,  high-order  solvers  for  geometrically  complex  configurations  containing  dielec¬ 
tric/magnetic  penetrable  scatterers,  open  and/or  closed  surfaces  as  well  as  perfect  and  lossy  conductors. 
These  are  configurations  of  fundamental  importance  in  diverse  fields,  with  application  to  (a)  Electro¬ 
magnetic  compatibility  (EMC),  (b)  Electromagnetic  interference  on  cavity-bound  electronics  (EMI),  (c) 
Evaluation  of  electromagnetic  response  of  dielectric/magnetic  coated  conductors,  and  (d)  Evaluation  of 
scattering  by  modern  metallic/nonmetallic  aircraft  structures — amongst  many  others. 


Figure  1:  Electromagnetic  scattering  by  a  perfectly  conducting  cubic  scatterer  (edge-length  =  2  under  plane- wave, 
k  =  20  normal  incidence.  Left:  ^-component  of  the  electric  field.  Center  and  Right:  x-  and  ^-components  of  the 
surface  current. 


2  Accomplishments 

Over  the  nine-month  duration  of  this  effort  we  developed  a  variety  of  electromagnetic  scattering  solvers 
whose  combined  use  enables  solution  of  a  wide  range  of  problems  in  the  field  of  electromagnetic  compat¬ 
ibility.  In  particular,  we  1)  Developed  surface  integral  equations  for  homogeneous  and  isotropic  dielectric 
bodies  whose  bounding  surfaces  can  contain  corners  and  edges ,  and  that  incorporate  regularizations  which 
give  rise  to  favorable  eigenvalue  distributions  and  small  numbers  of  GMRES  iterations;  2)  Implemented 
a  fast  high-order  solver  for  dielectric-body  integral  equations  introduced  per  point  1  above,  for  dielectric 


bodies  containing  edges  and  corners,  allowing  for  use  of  overlapping  and  non- overlapping  patches ;  3)  De¬ 
veloped  new  EM  scattering  solvers  for  open  metallic  surfaces  for  scatterers  containing  edges  and  corners; 
4)  Initiated  a  study  to  determine  the  domains  of  applicability  of  impedance  boundary  conditions  through 
comparisons  with  full  numerical  solutions  for  thin  volumetric  conductors;  and  5)  Produced  a  new  Green- 
function/Integral-Equation  methodology  for  solution  of  problems  involving  two  dimensional  periodicity 
in  three-dimensional  space;  see  e.g.  Figure  2  and  Section  2.3.  This  is  the  first  approach  ever  developed 
that  can  successfully  solve  bi-periodic  scattering  problems  in  three-dimensional  space  at  and  around  Wood 
anomaly  frequencies. 

Finally,  the  qualities  of  all  solvers  developed  were  demonstrated  via  compelling  applications  to  config¬ 
urations  relevant  to  the  field  of  electromagnetic  compatibility.  We  feel  we  have  successfully  addressed  the 
central  objectives  laid  down  on  our  Phase  I  proposal.  We  have  recently  received  an  invitation  to  submit  a 
Phase  II  proposal;  we  believe  the  results  of  our  work  over  the  Phase  I  effort,  which  are  described  in  what 
follows,  clearly  demonstrate  the  significant  potential  and  feasibility  of  the  proposed  methodologies. 

The  subsequent  Phase  II  work  will  proceed  to  1)  Enable  application  of  the  methodologies  developed 
during  the  Phase  I  to  multi-scale  geometries  containing  a  span  of  length-scales,  including  configurations 
such  as  full  aircraft,  attachments,  associated  electronics,  etc.;  and  2)  Produce  a  user  friendly  Graphical  User 
Interface  that  facilitates  application  of  the  software  infrastructure  to  relevant  engineering  configurations. 


Figure  2:  Scattering  by  an  infinite  bi-periodic  array  of  plates  with  periods  equal  to  ten  point  four  wavelengths: 
Pi  =  P2  =  10. 4A.  Plane  wave  incidence,  with  incidence  angle  0  =  45°.  Left:  21  x  21  section  of  the  infinite  array  and 
total  field  (incident  plus  scattered)  on  to  two  coordinate  planes.  Right:  Close-up  on  the  total  field  around  a  set  of 
three  plates. 


2.1  Well  conditioned  integral  equation  for  problems  of  scattering  by  dielectric  bodies 
and  thin  metallic  plates 

The  problem  of  scattering  by  dielectrics  is  governed  by  Maxwell’s  equations 

curl  curl  uJ  +  kj  uJ  =0  in  the  domains  ,  j  =  1,2 

for  the  electromagnetic  fields  uJ,  j  —  1,2  in  the  (bounded)  scatterer  D 2  and  the  surrounding  space  D1, 
respectively;  in  what  follows  we  assume  the  domains  D 1  and  D 2  are  separated  by  a  bounded  surface 
r.  The  material  properties  of  the  regions  DJ  are  characterized  by  two  wave-numbers  kj.  The  Maxwell 
equations  are  supplemented  with  (1)  transmission  boundary  conditions  on  the  material  interface  T  of  the 
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Figure  3:  Left.  General  composite,  consisting  of  a  number  of  penetrable  regions  containing  dielectric/magnetic 
materials  of  various  permittivities  and/or  permeabilities,  and  partially  coated  by  thin  metallic  surface  coatings  and/or 
containing  embedded  metallic  plates  or  bodies.  Right.  An  example  of  such  a  configuration:  Pendry’s  negative  index 
composite.  Manifold  important  applications  of  the  proposed  dielectric-metallic  electromagnetic  solvers  exist;  see 
points  (a)  through  (d)  above. 

form  y^u1  —  y^u2  =  — y^u mc  and  yjyii1  —  y)yU2  =  — y^umc  in  terms  of  the  Dirichlet  7^  =  n  x  •  and 
Neumann  traces  7^  =  n  x  curl  •  for  an  incident  field  umc,  and  (2)  radiation  conditions  at  infinity  on  u1. 
The  transmission  problem  in  dielectric  materials  with  bounded  interfaces  is  amenable  to  a  wider  variety  of 
boundary  integral  equation  formulations  than  is  the  corresponding  PEC  problem.  Among  these,  integral 
equations  of  the  second  kind  have  been  obtained  using  subtractive  cancellations  of  the  principal  parts  of 
hypersingular  operators  in  the  form  7/PV  —  (where  the  operators  Tj  correspond  to  the  wavenumbers 
kj).  However,  these  formulations  become  ill-conditioned  for  low  and  high-contrast  dielectrics  (see  also 
our  results  in  Figure  4)  and  thus  are  of  limited  use  for  a  range  of  applications  that  involve  such  material 
properties  (such  as,  e.g.,  narrow-band  dielectric  photonic  crystals). 

Since  the  availability  of  robust  integral  equation  formulations  throughout  the  range  of  dielectric  con¬ 
trasts  is  highly  desirable  for  treatment  of  the  types  of  problems  under  consideration,  we  seek  to  de¬ 
sign  such  formulations  on  the  basis  of  coercive  approximations  of  the  dielectric  admittance/capacity 
operators.  The  latter  operators  -reminiscent  of  the  Dirichlet-to-Neumann  acoustic  operators,  and  de¬ 
noted  hereby  as  7£T-are  defined  to  map  the  boundary  conditions  of  the  dielectric  transmission  prob¬ 
lems  for  the  field  u  =  (ux,u2),  that  is  yr(u)  =  (7^  —  yfjjyJv  ~  In)  (u)  to  the  Cauchy  data  of  u  on 
T  defined  as  7 c(u)  =  (y^u1,  y^u1;  y^u2,  y^u2),  and  thus  1Zt^t  —  7c-  Since  the  electromagnetic 
field  u  can  be  obtained  from  knowledge  of  its  traces  on  T  via  the  Stratton-Chu  formulas  in  the  form 
uJ  =  cj(jbu3  ,7wui)  =  (-l)3(MjtfDu3)  +  l/k?  £j{ y^iF)),  where  (Mj a)(z)  =  curl  fr  Gkj  (z-y)a(y)d<r(y) 
and  (<fya)(z)  =  curl  curl  fr  Gkj  (z  —  y)a(y )dcr(y),  it  follows  that  7t(Ci,  —  I  on  T.  The  main  idea  of 

our  approach  in  this  context  is  to  seek  suitable  approximations  of  the  admittance  operator  7 ZT  in  the  form 
(TZj ,  I—lZi)  — where  the  operator  approximates  in  a  certain  sense  the  operator  7^fy t  =  (7zn  7Jv) — and 

seek  for  electromagnetic  fields  uJ  ,  j  —  1,2  in  the  form 

u1  =  [Ci(nJ )](w),  u2  =  [C2(I  -  ^f)](w) 

for  an  unknown  tangential  current  distribution  w  =  (wi,  W2).  The  representation  presented  above  for  the 
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Figure  4:  Left:  number  of  GMRES  iterations  required  by  the  various  formulations  for  the  nearly  resonant 
electromagnetic  transmission  problem,  uo  =  50,  e\  =  1,  e^  =  2,  spherical  geometries-Diel-CFIE-RC  (red),  - 
Diel-CFIE-RPS  (blue),  -  Muller’s  CFIE  (black);  Right:  Number  of  iterations  required  by  each  formulation 
as  a  function  of  the  frequency  a;  =  1,1.1,. ..,50. 

fields  uJ  ,  j  —  1,2  leads  to  boundary  integral  equations  of  the  form 

[7T(Ci,C2)(^,/-^t)](w)  =  -(7tu-c,7>-c).  (1) 

Adopting  the  approximation  convention  that  1Z  ~  1Z  if  the  difference  TZ  —  TZ  is  a  compact  operator,  we  use 
the  following  operator  in  equation  (1) 

f,T  1  (  2fc2(|Ti  %-{)  i(kiT1  +  k2T2)  \ 

1  kj  +  qyihhih^  +  hTi)  2k\{fTiT2-{)  )  y) 

whose  derivation  is  based  on  the  pseudo-differential  calculus.  Taking  into  account  the  fact  that  7t(Ci,C2) 
is  the  exact  inverse  of  7 ZT ,  the  choice  of  IZj  in  equation  (2)  will  yield  regularized  integral  equations  (1) 
whose  spectra  will  accumulate  at  1.  Thus,  the  operators  TZ \  can  be  thought  of  as  regularizing  operators, 
and  the  ensuing  integral  equations  (1)  as  dielectric  CFIE-R  formulations.  Taking  a  cue  from  our  recent 
work  [1]  on  CFIE-R  formulations  for  the  PEC  case,  in  the  lossless  case  when  the  wavenumbers  kj  are  real 
we  will  complexify  them  in  the  form  kj  kj  +  iej,  6j  >  0  in  equation  (2)  in  order  to  obtain  coercive 
operators  TZj  which  thus  render  uniquely  solvable  integral  equations  (1). 

Results  of  our  ongoing  effort  suggest  that  significant  gains  can  result  from  use  of  our  proposed  regu¬ 
larized  formulations  (1).  Indeed,  in  Figure  4  we  compare  the  number  of  iterations  needed  for  a  GMRES 
residual  of  10-4  by  the  classical  Muller’s  formulations,  and  the  regularized  formulations  (1)  with  the  electric 
field  operators  replaced  by  their  principal  symbols  in  the  definition  of  the  regularizing  operators  TZ \ — the 
resulting  formulations  are  denoted  by  Diel-CFIE-RPS — and  the  complexified  wavenumbers — the  resulting 
formulations  are  denoted  by  Diel-CFIE-RC.  Clearly,  the  Diel-CFIE-RC  equations  consistently  produce 
significant  gains  in  numbers  of  iterations  throughout  the  frequency  spectrum  (Figure  4  right)  — about  one 
order  of  magnitude  reductions  compared  to  those  required  by  the  classical  Muller  formulations,  even  near 
resonances.  We  anticipate  that  commensurate  gains  will  be  produced  in  the  case  of  a  variety  of  dielec¬ 
tric  configurations  of  practical  interest  through  the  use  of  our  proposed  dielectric  CFIE-R  formulations. 
Furthermore,  our  methodologies  for  open  surfaces  and  for  singular  geometries  allows  for  a  seamless  incor¬ 
poration  of  the  dielectric  CFIE-R  formulations  in  the  case  where  the  interfaces  of  material  discontinuities 
exhibit  geometric  singularities  and  may  contain  lossy  or  lossless  metallic  components. 
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In  view  of  the  normal  and  tangential  singularities  of  surface  currents  for  open  surfaces,  we  seek  the 
scattered  electric  field  off  an  open  Perfect  Electrically  Conducting  (PEC)  surface  T  in  the  form 

E*(z)  =ik  J  Gk(\z  -  y\)Wl(y)ds(y)  +  V  J  Gk(\z  -  y\)divr(WI)ds(y),  (3) 

where  W  is  a  “weight”  matrix  and  where,  denoting  by  J  the  actual  surface  current  and  letting  J  =  W I 
(I  is  a  “regularized”  bounded  surface  current  and  W  carries  explicitly  the  current  singularity  at  the  edge), 
the  regularized  current  /  satisfies  the  EFIE  equation 


Twl  =  -n  x  Einc. 

Here  n  is  the  normal  to  T.  The  weighted  integral  operator  can  be  written  as 

Tw  I  —  <S\v  I  +  I 

where  operators  Sw  and  T>w  are  defined  as 

SWI  =  ikn  x  Gk(\x  -  y\)WI(y)ds(y) 


and 


T>w  I  —  — — curlp 
k 


J r  Gk(\x  -  y\)dWr(WI)ds(y) 


(4) 

(5) 

(6) 

(7) 


respectively.  Here,  for  a  scalar  field  F  defined  on  the  surface  T,  curlr  denotes  the  operator  curlrE  = 
(VF)  x  n. 

Let  r(u,v)  be  a  local  coordinate  chart  with  (u,v)  G  (0, 1)  x  (0, 1)  such  that  v  =  0  correspond  to  an 
open  edge.  In  this  case,  the  unknown  density  I  can  be  decomposed  as 


I  =  Iuru  +  Ivfv 


where  Iu  and  Iv  are  the  tangential  and  normal  components  of  the  surface  current.  Given  that  the  tangent 
and  normal  components  of  the  solution  have  singularity  of  0(1 /yfd)  and  0(y/d)  respectively,  where  d 
denotes  the  distance  to  the  edge,  the  choice  of  weight  matrix 


UJ 


1  -6cj2 

0  c 'jj 2 


(8) 


with  uo  ~  yfd,  and  6  —  ru  •  rv/ru  •  rn,  that  act  on  tangential  and  normal  components  ( IU,IV)T  of  I  as  a 
multiplication  by  the  matrix  renders  Iu  and  Iv  smooth. 

Before  we  discuss  the  numerical  method  that  is  used  of  to  solve  (4),  we  note  that  once  the  surface 
current  density  J  is  obtained  as  a  solution  of  (4),  the  electric  field  can  be  retrieved  using  equation  (3). 
Also,  if  we  define  the  electric  far  field  Eoo(x)  as 

Es(®)  =  t-(EooOk)  +  (9(M-1))  (9) 

\x\ 

then,  one  can  use  the  expression 

ju  r 

E00(x)  =  —xx  ^e-ik*y(WI(y)xx)ds(y)  (10) 

for  the  electric  far  field  computation. 
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Figure  5:  Electromagnetic  scattering  by  a  structure  made  out  of  an  array  of  5  x  5  thin  metallic  surfaces 
containing  corners  and  edges.  (As  discussed  in  the  text,  related  solutions  for  the  significantly  simpler  scalar 
case  were  available  previous  to  this  work.)  The  edges  are  treated  by  means  of  the  algorithm  described 
in  Section  2.2.1,  and  the  corners  are  modeled  by  extremely  sharp  rounded  edges,  as  shown  in  Figure  6. 
Structures  of  this  type  are  ubiquitous  amongst  the  configurations  and  devices  that  we  seek  to  treat  as  part 
of  the  present  effort;  compare  Figure  3. 

2.2  Electromagnetic  solvers  for  open  surfaces  with  corners 

These  methods  build  upon  previous  solvers,  which  were  described  as  part  of  our  STTR  proposal,  for  scalar 
scattering  by  open  surfaces  containing  corners.  The  present  solvers  provide  a  significant  generalization 
to  the  electromagnetic  case:  the  methods  and  implementations  for  solution  electromagnetic  scattering  by 
open  surfaces  containing  corners  is  the  result  of  work  associated  with  the  present  effort. 

The  presence  of  corners  in  a  scattering  geometry  results  in  solutions  that  become  unbounded  where 
this  blow-up  is  of  different  nature  from  that  at  open  edges.  It,  thus,  poses  a  significant  challenge  for  the 
numerical  scheme  described  above  in  term  of  achieving  high  accuracy.  In  view  of  this,  we  adopt  a  strategy 
where  we  solve  “nearby”  problems  on  smooth  domains  that  coincide  with  the  original  scatterer  except  in 
small  neighborhoods  near  the  corner.  The  success  of  the  strategy  hinges  on  its  ability  to  use  extremely  close 
approximations  of  the  domains  with  corners  and  use  of  discretization  that  give  rise  to  solution  with  limited 
computational  cost.  Here,  corners  have  been  smoothed  via  a  systematic  blending  of  two  arcs  on  either 
side  of  the  corner.  In  the  following  discussion,  we  describe  the  procedure  that  we  adopt  to  smoothly  round 
these  corners.  We  emphasize  that  the  procedure  is  completely  general  and  automatic:  it  can  be  applied  with 
ease  to  a  general  open  surface ,  with  curved  edges  and  associated  corners. 

2.2.1  Numerical  scheme 

In  this  section,  we  present  the  main  algorithmic  components  of  the  numerical  scheme  that  we  employ  in 
solving  (4),  which  in  turn  reduces  to  accurate  evaluations  of  the  weighted  integral  operator  7 w  and  hence 
S\y  and  Dw. 

We  note  that  the  accurate  evaluation  of  Sw  in  (6)  entails  obtaining  ( IU,IV)T  from  I,  which  is  quite 
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Figure  6:  A  unit  square  with  smoothly  rounded  corners.  The  geometry  is  divided  into  interior,  edge  and 
corner  patches  and  are  discretized  as  shown  above. 

straightforward.  An  application  of  the  VF-matrix  on  ( IU,IV)T  leads  to  integrals  of  the  form 

Sco[(p}( x)  =  f  -Gk(\x  -  y\)(j)(y)ds(y),  (11) 

Jr  u 

where  <j)  is  smooth,  that  can  then  be  approximated  through  a  specialized  numerical  integration  method. 
Toward  this  end,  we  employ  a  high  order  quadrature  that  relies  on  a  smooth  cut-off  function  r\x  supported 
in  a  small  neighborhood  of  the  target  point  x  that  is  also  identically  equal  to  one,  i.e.,  r)x  =  1,  in  the 
immediate  neighborhood  of  x,  to  split  (11)  as 

Su[4>](x)=  j  -Gk(\x  -  y\)yx(y)<f>{y)ds(y)  +  [  -Gk(\x  -  y\)(l  -  r]x{y))4>{y)ds(y).  (12) 

Jvx? o  w  Jr  ^ 

Clearly,  the  second  integral  in  (12)  can  be  accurately  approximated  using  a  “Weighted  Clenshaw-Curtis” 
quadrature  of  the  form 

rl  1 

/  ~rf(v)dv  «  y2wnf(vn),  (13) 
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where  the  quadrature  points  are  given  by 

t,”=Ki+cos((n+§)  £)) 

and  quadrature  weights,  wn,  by 


^ n 


1 

N 


m= 0 


4 

4m2  —  1 


cos 


run 

Jv 


where  the  primed  sum  denotes  that  the  first  term  (m  =  0)  is  halved.  The  first  integral  in  (12),  on  the 
other  hand,  is  evaluated  by  changing  to  polar  coordinates,  where  p-integrals  are  performed  using  a  scaled 
version  of  (13)  and  0-integral  is  evaluated  using  regular  Clenshaw-Curtis  quadrature  for  each  piecewise 
smooth  0-interval. 

In  the  evaluation  of  D\y  in  (7),  however,  care  should  be  taken  when  dealing  with  the  term  divr(VFI)  so 
that  the  singular  weight  uo  gets  proper  treatment.  To  this  end,  we  use  the  formula  for  the  surface  divergence 
in  coordinates.  For  this,  we  use  the  fact  that  for  a  given  tangential  field  X  =  Xuru  +  Xvrv  on  the  patch 
r(u,v ),  the  formula  for  the  tangential  divergence  reads: 

divrX  =  X  [du(^gXu)  +  dv{^gXv)]  (14) 

where  g  =  EG  —  F2  is  the  Riemannian  metric  tensor,  with  E  =  fu  ■  fu,  F  =  ru  ■  rv  and  G  =  rv  ■  fv.  A 
straightforward  application  of  (14)  thus  yields 


divr(WI)  = 


+ 


1 

L) 


2  g 


(dug(Iu  -  6oj2Iv)  +  dvg  co2Iv) 


duIu  -  du6u2Iv  -  9u2dur  +  u2dvIv  +  | dv(io2)Iv 


(15) 


assuming  that  uo  depends  only  on  v.  From  the  expression  in  (15),  it  follows  that  the  integral 

j  Gk(\x  -  y\)divr(WI)ds(y) 

in  (7)  has  the  same  form  as  in  (11)  and  thus  can  be  evaluated  using  the  integration  scheme  described 
above. 

The  last  remaining  element  of  our  numerical  algorithm  pertains  to  derivative  computations  that  arise 
in  (15)  as  well  as  in  the  surface  curl  differential  operator  in  (7).  As  these  expressions  only  involve  partial 
derivatives  of  smooth  functions,  one  can  use  Chebyshev  polynomials  as  spectrally  accurate  functional 
approximations,  which  can  then  be  used  for  finite  difference  approximation  of  the  derivatives  of  these 
functions.  One  can  also  directly  differentiate  the  approximating  Chebyshev  series  to  obtain  necessary 
derivatives.  In  this  case,  the  loss  of  accuracy  near  patch  edges  can  be  controlled  by  restricting  the  degree 
of  Chebyshev  polynomials  to  a  moderate  number. 


2.2.2  Preconditioned  Equation 

In  order  to  solve  better  conditioned  integral  equations,  we  precondition  equation  (4)  on  the  left  by  means 
of  the  regularizing  operator  7 ^  defined  by 

%K  =  ikn  x  J^Gk(\x  -  y\)WK(y)ds(y) 

-  ^curlr  J  Gk(\x  -  y\)u(y)divr(K)ds(y) 

—  iSu/K  A 
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(16) 


More  precisely,  we  solve 


(17) 


%,  o  Tw(l)  =  -%(n  x  Einc) 
where  we  express  the  operator  on  the  left-hand  side  as 

%o  o  Tw  —  Sw  °  Sw  +  Sw  °  T^w  +  °  Sw-  (18) 

We  note  that  the  operator  V ^  can  be  evaluated  in  a  straightforward  manner  using  differentiation  and 
integration  methods  described  in  the  previous  section. 


2.2.3  Numerical  Results 

Some  of  our  first  results  obtained  by  this  methodology  are  displayed  in  Figure  6,  which  shows  the  solution 
of  a  problem  of  electromagnetic  scattering  by  a  structure  made  out  of  thin  metallic  surfaces  containing 
corners  and  edges.  Structures  of  this  type  are  ubiquitous  amongst  the  configurations  and  devices  that  we 
seek  to  treat  as  part  of  the  present  effort;  compare  e.g.  Figure  3. 


2.3  Integral  equations  based  on  periodic  Green’s  functions 

A  variety  of  integral  equations  for  the  electromagnetic  transmission  and  reflection  problem  exist,  including 
those  arising  from  the  Stratton-Chu  representation  formulas,  that  express  the  scattered  fields  in  terms  of 
the  physical  surface  currents.  Integral  equation  formulations  of  the  dielectric  scattering  problems  are  posed 
in  terms  of  the  magnetic  and  electric  field  integral  operators  Kk  and  7&,  which  map  tangential  fields  a  into 
tangential  fields,  and  are  defined  by 

a)(x)  =  n(x)  x  j ^  VyGfc(x  -  y)  x  a(y)dcr(y),  (19) 

and 

(7fca) (x)  =  ikSka  +  ^7)fVa  =  ifcn(x)  x  J^Gk(x-  y)a(y)d<r(y) 

+  ^n(x)  x  Vx  J  Gfc(x  -  y)divra(y)d<7(y),  (20) 


where,  in  the  non-periodic  case,  Gk  is  the  outgoing  fundamental  solution  of  the  Helmholtz  equation  cor¬ 
responding  to  the  wavenumber  /c,  G^(x,  y)  =  G&(x  —  y)  =  4^|x-y| ;  the  hyper-singular  integral  in  the 
definition  of  K  should  be  interpreted  in  the  sense  of  Cauchy  principal  value.  We  denote  the  corresponding 
periodic  operators,  that  is,  the  integral  operators  arising  from  the  periodic  Green’s  functions  by  the  super¬ 
script  per.  Integral  equation  formulations  of  periodic  dielectric  scattering  problems  assume  as  unknowns 
the  (interior)  magnetic  and  electric  currents  j2  =  —  n  x  H2  and  m2  =  n  x  E2  and  takes  on  the  form  ICFIE 
on  T: 


(  (n2  +  /xi)(2//i)  li  +  1ic2  -  vy 

t  icjfjilSF  -  icofi2S2  +  i(ioe1)-1(Tr,Per  -  TD 


— iwei<Sf3r  +  iue2S2  +  1(TLPV,per  —  7^pv) 

(f‘2  +  ei)(2ei)_1I  +  e2eplC2  -  /Cfer 


X 


h 

m2 


nxff 
n  x  Ez 


(21) 


where  the  wavenumbers  ki  are  defined  as  ki  =  u^eiHi  for  i  =  1,2.  For  the  case  of  smooth  interfaces  T, 
the  integral  operators  /Q  and  7^PV’per  _  7^PV  are  compact  operators  from  HdS(T)  to  itself.  However,  the 
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operators  S{  are  not  compact  on  Nevertheless,  if  one  views  the  operator  on  the  left-hand  side 

of  equation  (21)  as  an  operator  from  HS(TM(T ))  x  HS(TM( r)),  where  HS(TM(T ))  denotes  the  classical 
Sobolev  space  of  tangent  vector  field  on  T,  then  all  the  operators  <%,  /Q,  and  7^PV  —  7^PV  are  compact. 

The  significant  computational  challenge  associated  with  the  numerical  solution  of  equations  (21)  is 
related  to  the  efficient  evaluation  of  the  periodic  Green’s  function  Gpkr.  There  is  an  extensive  literature 
devoted  to  the  issues  related  to  the  efficient  evaluation  of  periodic  Green’s  function.  In  what  follows 
we  present  a  method  of  solution  of  such  integral  equations  based  on  three  main  elements:  (1)  use  of 
windowing/ cutoff  functions  to  approximate  the  infinite  series  representations  of  Gpkr  by  truncated,  rapidly 
convergent  sums;  (2)  Taylor  expansions  of  the  terms  in  the  series  representation  of  Gpkr  corresponding 
to  ( m2d2  +  n2d 2)z  >  N  that  achieve  a  separation  of  the  contributions  of  the  integration  points  x'  from 
the  contributions  of  the  target  points  x;  (3)  high-order  quadrature  methods  based  on  overlapping  or 
non-overlapping  Chebyshev  patches  are  used  to  integrate  the  terms  in  the  relevant  integral  operators 
corresponding  to  the  contributions  in  Gpkr  that  arise  from  modes  such  that  ( m2d2  +  n  2 ^2) 2  <  and  (4) 
accelerated  evaluations  based  on  periodic  arrays  of  equivalent  sources  and  FFT  convolutions. 


Truncated  sums  using  windowing  functions.  The  very  slow/conditional  convergence  of  the  periodic 
Green’s  function  has  been  extensively  discussed  in  the  literature  and  several  methods  to  accelerate  its 
convergence,  notably  the  Ewald’s  method,  have  been  proposed.  We  propose  a  novel  idea  for  fast  evaluations 
of  periodic  Green’s  functions:  we  use  a  smooth  windowing  function  x  such  that  y(t)  =  l,t  <  1  and 
x(t)  =  0,t  >  2  and  we  approximate  Gpkr  in  the  following  manner: 


lmn  j  ^iandi  ^i/3md2 


Gfr(x>x')  ~  J2  x.  L 

(mdi)2  +  (nd2)2<4L2 

x  Gk(x  1  -  (x[  +  ndi),X2  -  [x2  +  rad2),x 3  -  £3)  =  G^er,L(x,  x7), 


(22) 


where  dmn  =  ((mdi)2  +  (mi2)2) 2  •  The  following  result  establishes  the  super-algebraic  convergence  of  Gper’J 


to  Gpker: 


Theorem  2.1  (Bruno,  Shipman,  Turc,  Venakides)  If  k  is  not  a  Wood  anomaly,  that  is  if  k2  ^  a ^  +  /32 
for  all  (m,  n)  G  Z  x  Z  where  am  —  a  +  fin  —  P  +  ,  then  for  all  x  ^  x'  and  all  integers  p  such  that 

2  <  p 

|Gr(x,x')  -  Gfr’L(x,x')|  <  CL\~p. 

Figure  7  demonstrates  the  excellent  accuracies  arising  from  use  of  Theorem  2.1. 


Separable  variables  representations  of  non-adjacent  interactions.  In  order  to  further  accelerate 
the  evaluation  of  Gper,L ,  we  derive  Taylor  series  expansions  of  quantities  G^{x  1  —  +  nd\),X2  —  (x2  + 

md2),X3  — X3)  for  |m|,  \n\  >  N  »  1  in  terms  of  the  small  parameter  {m2d\-\-n2d\)~^  and  expressions  that 
involve  separable  variables  x  =  (aq,  rr2,  £3)  and  x7  =  (x[,  x2,  x'3).  To  this  end,  we  introduce  the  notations 


rmn  =  {(x  +  mf  +  (y  +  nf  +  z2) 2 ,  r^n  =  (mz  +  nz) 2 , 


Urn'll,  — 


(m,  n ) 


and  express 


T  mn 


r^nni  1  T  a\OL\  +  n2o2  +  a2  +  a2  +  a2)  —  r^n 

^ mn  T  (*^j  ?/)  ’  ^ mn  ~f~  f  5 


1  +  —  (ql\cl\  +  o2a2)  +  O 


(23) 
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k  =  2 .5  Bloch  =  2*  (0  ,  0 .)  fcy,z)  =  (0.4,  0 .7,  0 .6) 


1.  1.5  2.  2.5 


Figure  7:  log10  of  the  periodic  Green  function  errors  as  a  function  of  the  log10  of  the  number  of  periods  L  considered 
in  the  smooth  truncation  of  the  periodic  Green  function,  according  to  Theorem  2.1. 


where 


x  y  z  2m  2  n 

Qj\  n  5  ^2  n  i  rj  5  ^1  n  5  ^2  rj  • 

pU  p  U  p  U  p  U  pU 


Using  the  fact  that  for  small  e 


(1  +  e)a  =  1  + 

where  5i  =  1/2  and  £  >  1  we  get  that 

OO  J 

/  =  rmn  ^(ai“i  +  a2«2  +  a?  +  a\  +  o|)£  -  xr^n^i0!  +  a2a2). 


£=1 


It  follows  then  that 


eikrmn  _  e =  _| -  ik f  —  ~k2  f 2  + 

2 


=  e^Le^(^2/)‘4n(l 


Furthermore,  we  have 


d  =  (rL)  Hl+fl1).  5  =  X!^^aiai+a2Q:2  +  al+a2  +  a3)€> 


where  71  =  —1/2  and  7^+1  =  —  ^777^  7^  and  hence 


e eik{x,y)-u°m. 


rO 


aikr ° 


Jk(x,y)-u\ 


0 

'mri 


(1  +  9  +  h  +  gh ) 

p,^)1 


1  + 


[1*1/2] 

i  W'  kjmn 


(j* 0  hi|—  j 
j_0  V  mn/ 


(24) 


(25) 


(26) 


(27) 


(28) 
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where  the  coefficients  Cm,n  are  of  the  form  nPlmP 2 
instance,  we  have 


777,77,  CllC  U1  LHC  1U1  111  ~T~q - sp  _|_p  and  can  be  computed  explicitly  for  all  m,n.  For 

\rmn) 


r>(i,o,o),o  _  ai  r»(o,i,o),o  _  a<2  r<(o,o,i),o  _  n 

^ mn  2  ’  2  5  ^ 

c%r,’a  =  =  5aia2 

C»rW  =  | (1  -a?/2),C£W  =  |(1  -  0^/2),  CfiW  =  |,C«*  =  |(1  -  W2). 


Consequently,  we  obtain  that 


—  (x[  +  ndi),X2  —  (x2  +  771^2)^3  —  £3)  = 

/ 


4ttgL 


0iAj(aj'i  — aJi  ,^2  ~x2  )'umn 


i+  v  x-x'H-D'ifriE 


where 


(i,i')^(0,0) 


dmn  —  {^TTld i)  +  ( 770?2 )  )2?  ^mn  — 


]i±C  .  .  \ 

2  c1+1,J  ' 


^ mn 


(mdi,  no?2) 

dm.n. 


and  we  used  the  classical  notation  x1  =  X ^-'2  x3  • 


(29) 


(30) 

(31) 


Fast  evaluations  of  the  boundary  layer  potentials  with  periodic  Green’s  functions.  Based  on 
the  ideas  developed  in  the  previous  sections,  we  approximate  the  evaluations  of  the  single  layer  potential 
in  the  following  manner: 


J  Gfer(x,X/)Mper(x/)^(x/) 


^  ^  ^  ^  dynn  ^  ^iandi  ei/3md2 


\m\,\n\<N 


L 


J  Gk{x  1  -  (x[  +  ndi),X2  -  (x'2  +  md2),x 3  -  x3 )/j,per(-x,)ds(xl) 


+  X  * 

iV<|m|,|n|<2L/di 


\  pikd-mn 

mn  \  eiandi  ei(3md2  _ _ e-ik{x\,x2)-Umn 


L 


An  dr, 


X  Xumn3^per){xi,X2:XzY  \  X 

|i|,|i'|=0  \  1  /  .  ~ 


li+i'l 

2  ^i+i  J 

V' 777, 77, 


(^n)li+i'|-i 


where  the  quantities  i7,  /iper)  are  defined  as 

I(umn,i,nper)  =  J  etk(xi’x>^'Umn  (x[,  x’2,  x’3 )Vper(x,)rfs(x/)- 


(32) 


(33) 


We  obtain  similar  expressions  for  all  of  the  other  necessary  boundary  layer  potentials  in  equations  (21). 
In  that  case  of  open  surfaces,  the  unknown  integral  density  fiper  is  singular  and  we  have  /iper(x)  = 
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where  /ire9  is  a  smooth  function  and  for  a  point  x  E  T  the  weight  cj(x)  =  (d(x,  dr))2.  Thus  we  pose  the 
corresponding  integral  equations  in  terms  of  the  smooth  density  jire9  via  the  ansatz  above.  Taking  into 
account  the  approximations  (32),  we  solve  the  following  integral  equations: 


E  * 

\m\,\n\<N 


giandi  ^ i(3md2 


+ 


X 


[  Gk{x  1  -  (X  +  ndi),x2  -  04  +  md2),x 3  -  3/3) 

Jr  ) 


ds(x ) 


E  x 

N<\m\,\n\<2L/di 


pikdmn 

giandi  ^ i/3md2  _ _ 

47rdmn 


ik{x1,X2)-Umn 


T 

E  7(  ^ mm  b  l^reg)  (x\ 5  X3)  (  1)^ 

iii;,it'i=o 

—  exp(ifcd  •  x),  x  G  T 


/n4+i,,J 

^mn 


(■ 4n)|i+i/| 


(34) 


where 


Jr  <^VX  J 


(35) 


and  the  expressions  umn ,  dmn  and  are  defined  in  the  previous  Section. 

In  order  to  obtain  high-order  discretizations  of  the  integral  operators  that  involve  periodic  Green’s 
functions  we  use  our  previously  developed  strategy  based  on  overlapping  patches  for  the  cases  when  the 
solutions  of  the  integral  equations  are  smooth  and  non-overlapping  Chebyshev  patches  in  the  cases  when 
the  solutions  undergo  rapid  transitions  in  the  vicinity  of  geometric  singularities.  Our  approach  calls  for 
(1)  the  evaluation  of  far-away  smooth  interactions  —  that  is,  contributions  to  the  integrals  related  to 
integration  points  x'  that  are  sufficiently  separated  by  the  target  points  x  and  (2)  the  evaluation  of  nearby 
singular  interactions  —  that  is,  contributions  to  the  integrals  related  to  integration  points  x'  that  are  close 
to  the  target  points  x.  In  case  (2),  the  terms  G^{x  1  —  (x[  +  ndi),X2  —  ( x'2  +  —  x'3 )  are  singular 

only  for  x  =  x'  and  n  —  m  —  0 — that  is  the  free-space  Green’s  function,  and  the  singularity  is  resolved 
with  high-order  by  polar  changes  of  variables  and  interpolation  techniques.  The  far-away  interactions 
(1)  in  the  case  of  the  free-space  Green’s  function  Gk  are  accelerated  by  use  of  equivalent  sources  placed 
on  Cartesian  grids  and  3D  sparse  FFT  convolutions.  These  techniques  extend  to  the  case  of  evaluating 
far-away  interactions  via  G^per  which  we  describe  next. 

The  first  step  of  this  approach  consists  of  partitioning  a  cube  C  of  size  A  circumscribing  the  scatterer 
into  L3  identical  cubic  cells  C{  of  size  adjusted  (in  the  sense  of  small  enough)  so  that  they  do  not  admit 
either  inner  acoustic  resonances  -  eigenfunctions  of  the  Laplacian  with  Dirichlet  boundary  conditions.  The 
main  idea  of  the  acceleration  algorithm  is  to  seek  to  substitute  the  surface  “true”  sources  which  correspond 
to  the  discretization  points  contained  in  a  certain  cube  C{  by  acoustic  periodic  “equivalent  sources”  on  the 
faces  of  Ci  in  a  manner  such  that  the  acoustic  fields  generated  by  the  ^-equivalent  sources  approximate 
to  high  order  accuracy  the  fields  produced  by  the  true  C{  sources  at  all  points  in  space  non- adjacent  to 
Ci.  The  precise  concept  of  adjacency  [2]  results  from  requiring  that  the  approximation  corresponding  to 
a  given  cell  C{  be  valid  within  exponentially  small  errors  outside  the  concentric  cube  S{  of  triple  size.  At 
the  heart  of  this  method  lies  the  use  of  equivalent  sources  which  consist  of  acoustic  monopoles  and  dipoles 
placed  on  three  independent  sets  n^c,  each  one  parallel  to  xi  =  0.  For  a  fixed  value  l  —  1,  2,  3,  we  associate 
to  an  acoustic  field  u  and  each  cell  ^-equivalent  sources,  acoustic  monopoles  G^lPer(x  —  x[  j)  and 

dipoles  dG^per(x  —  x\  -)/dxi  placed  at  points  x\  -^  l  —  1,  •  •  •  ,  Mequiv  contained  within  certain  subsets 

which  lie  within  the  union  of  two  circular  domains  concentric  with  and  circumscribing  the  faces  of  q, 
their  radius  chosen  according  to  the  prescriptions  in  [2].  The  fields  ^>ciJrue  radiated  by  the  C{- true  sources 
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Scatterer 

Unknowns 

L 

e 

^ conv 

It /Time 

Sphere 

6  x  16  x  16 

20 

8.0  x  10-3 

1.8  x  10"2 

7/3ml2sec 

Sphere 

6  x  16  x  16 

30 

1.4  x  10“4 

1.8  x  10“4 

7/7m3sec 

Sphere 

6  x  16  x  16 

50 

2.8  x  10-5 

7/19m33sec 

Sphere(T  Ac) 

12  x  16  x  16 

20 

2.1  x  10-3 

8.2  x  10"4 

19/5m38sec 

Sphere(T  Ac) 

12  x  16  x  16 

30 

8.4  x  10“4 

2.4  x  10“4 

19/9m29sec 

Sphere(T  Ac) 

12  x  32  x  32 

50 

4.9  x  10-5 

19/26ml0sec 

Sphere(T  Diel) 

24  x  16  x  16 

20 

7.2  x  10-4 

9.8  x  10"5 

103/28m23sec 

Sphere(T  Diel) 

24  x  16  x  16 

30 

6.3  x  10-4 

8.0  x  10“6 

103/31m50sec 

Sphere(T  Diel) 

24  x  32  x  32 

50 

3.4  x  10-5 

103/44m8sec 

Cube 

6  x  16  x  16 

20 

1.6  x  10-3 

6.4  x  10"3 

34/7m57sec 

Cube 

6  x  16  x  16 

30 

8.6  x  10“5 

2.3  x  10“3 

34/llm42sec 

Cube 

6  x  16  x  16 

50 

5.5  x  10-5 

34/53ml3sec 

Disc 

5  x  16  x  16 

20 

1.3  x  10-3 

3.7  x  10"3 

10/4m37sec 

Disc 

5  x  16  x  16 

30 

7.1  x  10“5 

4.2  x  10“4 

10/8ml8sec 

Table  1:  Convergence  of  the  solvers  using  G 


L,per 

k 


even  without  use  of  Taylor  expansions . 


are  approximated  themselves  by  fields  ^pc^eq  radiated  by  the  Q-equi  valent  sources 

eq 


—  Me(lu'i'v 


'(*>=  E 

3= 1 


Jm)l  nL,per 

kj  G 


(x>  xi,o + if 


dGLk’per(x,x. 


dxi 


(36) 


The  parameters  Mequiv  and  the  unknown  monopole  and  dipole  intensities  in  (36)  are  chosen  so  that  the 
truncated  spherical  wave  expansions  of  order  nt  for  ^c^true  and  ipCi’eq  differ  in  no  more  than  0(e)  outside 
Si.  Based  on  considerations  on  spherical  harmonics,  it  was  required  in  [2]  that  Mequiv  >  n2  equivalent 
sources  are  used  for  each  acoustic  component  and  the  the  intensities  are  chosen  such  that  to  minimize  in 
the  mean-square  norm  the  differences  ('0Q,eg,(x)  —  'ipCi,true (x.))  as  x  varies  over  a  number  nco11  collocation 
points  on  dSi.  Hence,  the  intensities  in  (36)  are  obtained  in  practice  as  the  least-squares  solution  of  three 
overdetermined  linear  systems  =  b  where  A  are  nco11  x  Mequiv  matrices.  This  strategy  leads  to  a 
computational  cost  of  0(4L2A4/3  log  A)  to  evaluate  the  boundary  layer  potentials  involving  GGper,  where 
A  is  the  number  of  discretization  points. 

Finally,  the  terms  in  equations  (33)  and  (34)  that  involve  separate  contributions  of  the  target  point 
x  from  the  integration  point  x7  are  computed  in  the  following  manner:  (a)  for  all  (m,  n)  such  that  A  < 
|m|,|n|  <  2 L  the  coefficients  Cmn  are  precomputed  for  all  multi-indices  i  such  that  |i|  <  2 T  and  all 
indices  j  such  that  o  <  3  <  T  following  the  Taylor  series  algebra  manipulations  described  in  Section  3.1 
and  (b)  for  all  multi-indices  i  per  part  (a)  the  integrals  I(uj,  i,  ixreg)  are  computed  for  a  number  J  of 
directions  j  uniformly  distributed  on  the  two-dimensional  unit  circle  S1  and  the  quantities  I{umn,  i,  tireg) 
for  all  (m,  n)  such  that  A  <  |m|,  \n\  <  2 L,  in  turn,  are  obtained  from  Fourier  interpolation  from  the  values 
I{uj ,  i,  iireg ),  j  =  1, . . .  J —  the  number  J  is  typically  much  smaller  than  4(2L  — A)2.  Thus,  we  must  choose 
three  parameters  A,  L,  and  T  in  equations  (34)  so  that  we  get  high-order  accuracy  in  small  computational 
times. 


2.4  Numerical  results 

In  this  section  we  present  several  numerical  results  concerning  scattering  off  of  periodic  arrays  of  scatterers. 
In  order  to  assess  the  accuracy  of  our  solver,  we  present  the  energy  balance,  that  is  the  energy  of  the  incident 
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wave  equals  the  energy  of  the  reflected  wave  plus  the  energy  of  the  transmitted  wave.  The  energy  balance 
can  be  expressed  mathematically  in  terms  of  the  Rayleigh  coefficients  (also  called  “efficiencies”),  and  it 
takes  on  the  following  form: 


^  7rc,ra|^n,m  +  ^ra’ral  —  7*  (37) 

(n,m)£U  (n,m)eU 


in  the  acoustic  case  and 


7n,ra|B^m|  +  £  7n,m|Bn,m  +  ^n’ml  —  7-  (38) 

(n,m)Gt/  (n,m)EU 


in  the  electromagnetic  case.  We  will  present  in  our  numerical  tests  the  energy  balance  which  is  defined  as 


e  — 


^2(n,m)eU  7n,m|^n,m|2  +  J2(n,m)eU  7n,m|^n5m  +  ^n,m\2  7 


7 


(39) 


In  the  case  of  transmission  problems,  the  transmission  quantity  is  a  measure  of  the  fraction  of  the  energy 
of  the  incident  wave  which  is  transmitted  by  the  periodic  scatterer  and  it  is  defined  as: 


0,0  .2 


T  = 


7 


(40) 


We  present  numerical  results  for  three  geometries:  spheres,  cubes,  and  discs.  The  first  case  is  treated 
with  the  overlapping  method  and  the  other  two  with  the  non-overlapping  method.  We  present  the  energy 
balance  e  as  well  as  the  errors,  evaluated  by  means  of  a  resolution  study,  in  the  coefficient  £>q~0  in  the 
sound-soft  case  and  in  T  in  the  transmission  case-we  denote  these  errors  by  econv.  All  linear  systems  were 
solved  using  GMRES  to  a  residual  of  10-4. 

We  present  in  Table  1  the  high-order  nature  of  our  solver  as  a  function  of  the  parameter  L  in  the 
definition  of  G^per .  In  all  of  the  sound  soft  cases  the  periods  are  d\  —  —  4,  we  consider  normal 

incidence  ^  =  <j>  =  0,  the  wavenumber  k  =  0.75.  In  the  case  of  the  transmission  problem,  we  take  cj  =  1, 
=  1,  62  =  40  and  =  7r/6,  <j>  =  0.  With  regards  to  computational  times,  for  comparison  we  note  that 
the  implementation  one  of  the  most  advanced  techniques  for  evaluation  of  periodic  Green’s  functions  [3] 
(which  is  based  on  Kummer  transforms,  either  spatial  or  spectral  representations,  supplemented  by  Shanks 
transforms)  is  reported  to  take  several  milliseconds  per  evaluation  point  [4].  Thus,  for  a  discretization 
6  x  16  x  16  there  are  about  1.8  x  106  evaluations  of  periodic  Green’s  functions  which  will  require  at  least 
1.8  x  103  seconds  (30  minutes)  to  evaluate  just  one  matrix  vector  product.  In  contrast,  as  it  can  be  seen  in 
Tables  1-4,  our  solvers  require  about  82  sec  for  a  matrix  vector  product  to  obtain  results  with  four  digits 
of  accuracy. 

In  the  next  three  tables,  Tables  2-4  we  present  results  for  the  cases  when  Taylor  expansions  are  used. 
The  results  concern  the  sound-soft  case,  normal  incidence,  d\  —  g?2  =  4.  We  conclude  that  reductions  about 
one  order  in  magnitude  in  computational  times  can  be  garnered  from  use  of  Taylor  expansions,  while  the 
results  are  still  accurate  with  two  digits. 

We  present  in  Figure  8  the  acoustic  transmission  coefficient  as  a  function  of  the  frequency  cj  in  the  case 
of  scattering  from  a  periodic  array  of  spheres  with  62  =  40  and  two  values  of  the  azimuthal  angle  —  0 
(left)  and  =  7r/6  (center)  and  a  periodic  array  of  smoothed  cubes  (d=0.1)  (right)  with  normal  incidence 
and  62  =  40.  We  note  that  sharp  transitions  in  the  transmission  curves  correspond  to  resonant  behavior 
—  around  uo  —  0.7  for  normal  incidence  and  uo  —  0.5  for  the  case  of  oblique  incidence  in  the  case  of  spheres 
and  for  cj  =  0.79  in  the  case  of  smoothed  cubes. 
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k 

Unknowns 

N 

L 

T 

e 

^conv 

It /Time 

1.75 

1.75 

6  x  16  x  16 

6  x  16  x  16 

20 

14 

20 

20 

2 

2.8  x  10-2 

5.8  x  10“3 

3.5  x  10-2 

9/3ml5sec 

10/lm20sec 

3.75 

3.75 

6  x  16  x  16 

6  x  16  x  16 

80 

20 

80 

40 

2 

1.9  x  10“3 
1.4  x  10-2 

1.3  x  10-2 

18/51m4sec 

18/3m41sec 

4.75 

4.75 

6  x  16  x  16 

6  x  16  x  16 

80 

28 

80 

80 

2 

1.9  x  10-3 

5.9  x  10“3 

2.0  x  10“3 

29/54m51sec 

29/9m59sec 

Table  2:  Scattering  by  a  by-periodic  array  of  spheres;  e  and  econv  denote  the  energy  balance  error  and 
the  error  in  reflection/transmission  evaluated  by  means  of  a  convergence  study.  Note  that  use  of  a  Taylor 
expansion  (indicated  in  the  table  by  displaying  the  Taylor  truncation  order  T  —  2)  significantly  accelerates 
the  solution  process. 


k 

Unknowns 

N 

L 

T 

e 

£ conv 

It /Time 

1.75 

1.75 

6  x  16  x  16 

6  x  16  x  16 

20 

14 

20 

20 

2 

9.0  x  10“3 
2.0  x  10-2 

4.2  x  10-2 

50/10m4sec 

50/8m42sec 

3.75 

3.75 

6  x  16  x  16 

6  x  16  x  16 

80 

20 

80 

40 

2 

1.3  x  10-5 
1.9  x  10“3 

1.9  x  10“2 

70/58m42sec 

70/14m21sec 

4.75 

4.75 

6  x  16  x  16 

6  x  16  x  16 

80 

28 

80 

80 

2 

5.7  x  10“3 
9.0  x  10-3 

2.9  x  10-2 

80/60m8sec 

80/27m59sec 

Table  3:  Scattering  by  a  by-periodic  array  of  cubes;  otherwise  same  as  figure  2. 


k 

Unknowns 

N 

L 

T 

e 

^ conv 

It /Time 

1.75 

1.75 

6  x  16  x  16 

6  x  16  x  16 

20 

14 

20 

20 

2 

4.0  x  10“2 
1.2  x  10-2 

4.1  x  10-2 

ll/4m47sec 

ll/3m8sec 

3.75 

3.75 

6  x  16  x  16 

6  x  16  x  16 

80 

20 

80 

40 

2 

1.4  x  10-5 
2.0  x  10“2 

2.4  x  10“2 

17/50m56sec 

17/5m54sec 

4.75 

4.75 

6  x  16  x  16 

6  x  16  x  16 

80 

28 

80 

80 

2 

6.5  x  10“3 
2.3  x  10-2 

4.4  x  10-2 

20/51m8sec 

29/10m47sec 

Table  4:  Scattering  by  a  by-periodic  array  of  plates  (see  Figure  2);  otherwise  same  as  figure  2. 


Figure  8:  Transmission  coefficients  as  a  function  of  the  frequency  u  in  the  case  of  scattering  from  a  periodic 
array  of  [spheres  with  62  =  40  and  ijj  =  0  (left)  and  ^  =  7t/6  (center)  and  from  a  periodic  array  of  cubes 
(right). 
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To  conclude  this  report  we  return  to  its  first  image,  Figure  1,  which  displays  the  dielectric  transmission 
coefficient  as  a  function  of  the  frequency  u  in  the  case  of  scattering  from  a  periodic  array  of  spheres  with 
62  =  40  and  normal  incidence  around  the  resonant  frequency  u  —  0.69,  and  the  fields  inside  and  around  a 
spherical  element  in  a  periodic  array.  Experiments  such  as  these  demonstrate  the  ability  of  the  new  solvers 
to  produce  rapidly,  efficiently  and  accurately  important  physical  observables  such  as  total  reflection  and 
total  transmission  for  challenging  three  dimensional  photonic  crystals.  In  view  of  the  computational  times 
reported  above  in  this  text,  we  believe  these  solvers  are  orders  of  magnitude  faster,  for  a  given  accuracy, 
than  all  other  methods  available  at  present  and,  in  fact,  will  enable  solution  of  previously  intractable 
problems. 
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